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FAST APPROACH FOR CALCULATING FILM THICKNESSES AND PRESSURES 
IN ELASTOHYORODYNAMICALLY LUBRICATED CONTACTS AT HIGH LOADS 

Luc G. Houpert* and Bernard J. Hamrock 
National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


SUMMARY 

The film thicknesses and pressures In elastohydrodynamlcal ly lubricated 
contacts have been calculated for a line contact by using an Improved version 
of Okamura's approach. The new approach allows for lubricant compressibility, 
the use of Roelands' viscosity, a general mesh (nonconstant step), and accurate 
calculations of the elastic deformations. The new approach Is described, and 
the effects on film thickness, pressure, and pressure spike of each of the 
Improvements are discussed. Successful runs have been obtained at high pres- 
sure (to 4.8 GPa) with low CPU times. 
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SYMBOLS 

weighting factors used to define dP/dX at node 1 
half Hertzian length, R ^8W/w, m 
weighting factors used to Integrate P 

Influence coefficient used to calculate elastic deformation at 
node 1 due to Pj 

v oung's modulus of surface, Pa 


equivalent Young's modulus, 


1 

E' 



materials parameter, aE' 
dimensionless film thickness, hR/b 2 * *h/8RW 
dimensionless film thickness where dP/dX * 0 
dimensionless constant used In calculation of H 
maximum deviation from flat shape of Hertzian contact 
film thickness, m 


♦Presently at S.K.F. Engineering and Research Center, 3430AB Nleuwegeln, 
The Netherlands; NRC-NASA Research Associate during academic year 1984-85. 



nodes 


U 

N number of nodes used In linear system 

N max maximum number of nodes used In mesh 
P dimensionless pressure, p/p^, Pa 

P s dimensionless pressure at spike 

p pressure 

PH maximum Hertzian pressure, E'b/4R » E' ^U/2t, Pa 
R equivalent radius of contact, m 

r radius of surface, m 

U dimensionless speed parameter, nou/E'R 

u tangential speed of surface, m/s 

u average entrainment rolling speed, (u a ♦ Ub)/2, m/s 
M dimensionless load parameter, w/E'R 

w applied load per unit length, N/m 

X dimensionless abscissa, x/b 

X er>( j outlet meniscus distance 

*max maximum value of X In mesh 

x m1n minimum value of X In mesh, X-| 

ax constant step, (X^ - X-j)/(N - 1), used In uniform mesh 
x abscissa along rolling direction, m 

Z Roelands' parameter 

a plezovlscous coefficient, m 2 /N 

maximum dimensionless Hertzian deformation 
4^ elastic deformation at node 1, m 

4^ dimensionless elastic deformation, R4/b 2 

no viscosity at operating temperature and ambient pressure, N s/m 2 

v Poisson's ratio 
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p relative density 

p e relative density where H * H e 

Subscripts: 


a 

surface 

a 

B 

Barus 


b 

surface 

b 

H 

Hertz 


1 

at node 

number 1 

J 

at node 

number J 

max 

maximum 

value 

min 

minimum 

value 

R 

Roelands 


INI RODUCTION 

Elastohydrodynamlc lubrication Is a form of fluid-film lubrication where 
elastic deformations of the lubricated surfaces become significant. It Is usu- 
ally associated with highly stressed machine elements such as rolling-element 
bearings and gei.-s. The first attempt to analyze elastohydrodynamlcally lubri- 
cated contacts was made by Grubln and Vinogradova (ref. 1), who managed to 
Incorporate both the effects of elastic deformation and the viscosity-pressure 
characteristics of the lubricant In an Inlet analysis of hydrodynamic lubri- 
cation of nonconformal contacts. Solutions to the coupled Reynolds and elas- 
ticity equations have proved to be difficult. An attempt to review the methods 
of obtaining numerical solutions to the elastohydrodynamlc lubrication problem 
Is presented by Hamrrck and Tripp (ref. 2). Highlights of four main approaches 
are presented, namely the direct method, the Inverse method, the quasi-inverse 
method, and the system approach method. The advantages and disadvantages of 
each method are discussed. As described by Hamrock and Jacobson (ref. 3), the 
direct method used a large amount of CPU time and the approach failed at high 
loads. (Maximum dimensionless load W was of the order of 3xl0” 5 ). 

A recently developed approach not covered In Hamrock and Tripp's review 
(ref. 2) Is that of Okamura (ref. 4). Okamura uses a Newton-Raphson method to 
solve the elastohydrodynamlc lubrication problem In a system form. Iterations 
on the system approach are still required, as was the case In the direct 
method; however, much less computer time Is used. Furthermore, the major 
advantage of the method Is that solutions can be obtained for higher loads than 
by the direct method. (Maximum dimensionless load W was of the order of 
8xl0~ 5 .) Nevertheless, these loads are much lower than those In rolling- 
element bearings, where W reaches the order of lxlO -3 . At high loads the 
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elastic deformations are several orders of magnitude larger than the film 
thickness. Any Inaccuracy In the pressure will cause (via elasticity calcula- 
tions) drastic changes In the film thickness that will In turn (via the 
Reynolds equation) cause large pressure fluctuations and numerical Instabil- 
ities, especially when the viscosity Is large, as Is the case In normal EHL 
contacts. At high loads the viscosity of the fluid can vary by 10 orders of 
magnitude within the conjunction. 

In the current study It was felt that the only way of solving the EHL 
problem at high loads was to accurately calculate the elastic deformations and 
the pressure gradient dP/dX, especially In the Inlet region and near the pres- 
sure spike, where dP/dX Is large. To calculate the pressure gradient, a 
general mesh (l.e., nonconstant step) Is needed (a fine step when dP/dX Is 
large). Okamura's method (ref. 4) was therefore extended to Include these 
Improvements . 

For a comparable change In pressure the relative density change Is con- 
siderably smaller than the viscosity change. The maximum density Increase Is 
about 35 percent. However, although the density effect Is much smaller than 
the viscosity effect, the very high pressures In EHL films are such that the 
liquid should no longer be considered Incompressible. Furthermore, high- 
pressure viscosity measurements have shown that Barus 1 expression cannot be 
used at high loads. A more realistic pressure-viscosity relationship as pro- 
posed by Roelands (ref. 5), should be used. Although Roelands expression Is 
valid for pressures not exceeding 1 GPa, It can be used In EHL calculations at 
larger maximum pressures because the film thickness Is mainly dependent on the 
Inlet of the contact, where pressures are small. Therefore, lubricant compres- 
sibility and Roelands* viscosity were also Introduced In the present approach. 
Different means (Including the authors' new approach) of calculating the 
elastic deformations are first compared. After having described the new mathe- 
matical formulation, the effects on film thickness, pressure, and pressure 
spike of the mesh and number of nodes are studied In detail. Results are shown 
for low, Intermediate, high, and very high loads. Finally results obtained 
using Barus* and Roelands* viscosity as well as the compressible and Incom- 
pressible relationships are compared. 

The authors wish to thank Dr. J. Tripp of the NASA Lewis Research Center 
and Mr. T. Lubrecht of the University of Twente, The Netherlands, who took 
part In many helpful discussions on the subject. 


CALCULATION OF ELASTIC DEFORMATIONS 

The elastic deformation 6 at any point x on the surface Is expressed 
as 
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In (x - x') 2 dx* 


( 1 ) 


where p Is the pressure function of x_|_, with x* varying between x m ^ n 
and x en( j. Letting X - x/b, P » p/Ph. 4 » R4/b 2 where 
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( 2 ) 


1/R - 1/R a ♦ 1/R b , R/b -yt/8W, U - w/E'R, and p H » E'b/4R gives 


lm ~U 


/* X end 2 

■) 

/* X end 

J P In (X - X') dX* ♦ In 

1 PdX ' 

, X m1n 


A m1n 


But w « y*p dx Implies /p dX 1 * «/2. Making use of this, equation (2) 
becomes 




end 


P In (X - X') 2 dX' 


min 


- ! " (>’ ?) 


(3) 


The constant term on the right side of equation (3) represents In general 80 
to 90 percent of the total deformation. This Is a very useful separation In 
that It puts the remaining pressure-dependent deformation of the same order as 
the film thickness at moderate loads. Using Integration by parts on the 
Integral gives 



where 

I * y* In (X - X') 2 dX ‘ - -(X - X') [In (X - X') 2 - 2] (5) 

Since the pressure Is zero at X m ^ n and X en( j, we can write that 

i * - £ f Cnd jjJ.tX - X 1 ) [In (X - X 1 ) 2 - 2] dX* - \ In ^R 2 (6) 

X m1n ' ' 

The Integral In equation (6) can be calculated analytically by assuming that 
the pressure Is described by a polynomial of second degree In the Interval 
[Xl-1* x Hll- The details of the calculations are given In appendix A. The 
deformation at node 1 Is calculated as a function of the pressure Pj and 
the Influence coefficients D^j: 

Vi- i ’* (•*!“) (7) 

The results obtained will be compared In table I with the ones obtained by 
Hamrock and Jacobson (ref. 2) and Okamura (ref. 4), who used simpler approaches. 
Hamrock and Jacobson assumed the pressure to be constant In the Interval 
[Xj - AX/2, Xj ♦ AX/2] and used an analytical expression for f he Integral of 
log (Xi - X'). Okamura did not use any analytical solutions and assumed 
simply that 
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(8) 



The three approaches can be compared by assuming a Hertzian pressure. Between 
abscissa -1 and *1 , the Hertzian shape should be flat, leading to AH * 0 where 
AH Is defined as 


AH 



( 9 ) 


and 6 max Is the maximum deformation and can be compared with the Hertzian 
analytical value Jh: 






♦ \ In (2) 


0.25 


( 10 ) 


The largest value AH^x/Hg of AH/H e Is found at abscissa X » -1 
and *1 (because of thejslope^ discontinuity) and Is shown on table 1 with the 
corresponding value of imax^^H * 1* The central film thickness H e has 
been chosen to be equal to 0.5; X m ^ n and X max define the first and last 
values of X; N^x Is the number of nodes. 


TABLE I. - THREE WAYS OF CALCULATING ELASTIC DEFORMATIONS 


[OK: O'.amura (ref. 4); HJ: Hamrock and Jacobson (ref. 3).] 


tynax 

Xfnin 

X max 

( 4 max/ 4 H “ 1 1 

4 *1nax^e 

OK 

HJ 

New approach 

OK 

HJ 

New approach 

i>l 

-1.0 

1.0 

-1.9x10-3 

-2.7xl0 -3 

8.6x10-7 

9.4x10“ 3 

-5.5xir 3 

-2.7xlO" 3 

bl 

-3.b 

1.4 

-4.9xl0 -3 

-l.ixicr 2 

l.OxlO" 5 

3.4xlO" Z 

-1.6xlO -Z 

-8.6xl0“ 3 

ibl 

-1.0 

1.0 

-6.3x10-!* 

-5.1xl0 -4 

5.1x10-8 

4 .OxlO -3 

-1.4xl0“ 3 

-6.4xl0" 4 

151 

-3.6 

1.4 

-1.6xl0“ 3 

-2.0xl0“ 3 

5.4x10-' 

1.2xl0" z 

-4.4xl0" 3 

-2.1xir 3 

301 

-1.0 

1.0 

-3.ixicr 4 

-1.8xl0“ 4 

8.1x10-9 

2.2x10“ 3 

-5.6xl0" 4 

-2.5xlO“ 4 

301 

-3.6 

1.4 

-7.8xl0“ 4 

-7.ixicr 4 

9.0xl0" 8 

6. OxlO" 3 

-1.8xir 3 

-8.6xl0r 4 

6b 1 




2.2xl0" 6 



-2.6xl0“ 4 


-3 . o 

1 .H 





51 


*1.0 



6.5xl0- 7 



-2.7xlCT 4 

#U 







a Nonuniform. 


The results of table I show that for a given mesh the best accuracy In 
the calculation of A and AH Is obtained by using the present approach, 
lhe value of 4 max from Hamrock and Jacobson (ref. 3) Is In some cases less 
accurate than that from Okamura (ref. 4) because they did not use the constant 
expressed In equation (3). But Is not really significant since any 

Inaccuracy In Its calculation can be compensated for by Hq In the film 
thickness expression. The Important parameter of table I Is AH max since 
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It Is a measure of the flatness of the film shape. This aspect Is extremely 
Important at high loads, where the elastic deformations are two or three orders 
of magnitude larger than the film thickness. 

Also shown In table I are the results obtained with a uniform mesh of 661 
nodes by Hamrock and Jacobson (ref. 3). Very small values of aH max were 
calculated and cannot be reproduced with the new approach because of storage 
problems with the matrix and because of the large system of equations 

that would have to be solved with such a mesh. But using a nonuniform mesh 
with a fine grid near X equal to -1 and +1 (where 10 sublrtervals are used In 
the Intervals [-1 ,-0.96] and [0.96,1]), similar results are found with a very 
small value of N^x (N^x ■ 51), as Indicated In table I. The latter case 
Illustrates very well the power of the new approach. 


CALCULATION OF dP/dX 


The term dP/dX Is used to calculate the Influence coefficient O^j 
and Is used In the Reynolds equation described later. For the calculation of 
D^j, (dP/dX)^ at node 1 Is calculated by using three nodes. 

(dx) 1 “ a 1 ,1-1 P 1-l * a 1,1 P 1 * a 1,1+l P 1+l (1 


where 


X - X 

*1 A 1+l 

•1.1-1 " (x 1 _ 1 - x 1 )(x 11 - x ul ) 

2X i - x i+i ~ x i-i 

* 1.1 * (X^ - X 1 _ 1 )(X 1 - x ul ) 

X 1 - x 1-l 

1 * 1*1 “ < X W 1 - X 1 - 1 >< X 1+1 - X 1 > 


( 12 ) 

(13) 

(14) 


At the first and last nodes, X 1 and X u , dP/dX Is also defined by 

max 

using three nodes (X 1 to X 3 and X N? to X^ ). 

max 

For the calculation of dP/dX In the Reynolds equation, two nodes may 
often be used If the numerical convergence Is difficult to obtain with the 
three-node formula. For the two-node formula the weighting factors are 


'1,1+1 


X _ X ’ 
* 1+1 A 1 


a 1 ,1 " " a 1,1+l ; 


‘ 1 , 1-1 * 


( 15 ) 


Minor corrections are also applied on the last values of a n j to respect 
the boundary conditions at X * X en( j and are described In appendix B. 
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FORMULATION OF MATHEMATICAL PROBLEM 


The basic equation to solve In dimensionless form at each node 1 Is the 
dimensionless Reynolds equations. 


with 




0.6xl0" 9 P H P 1 

^1 - 1 ♦ Iq 

1 f 1.7x10 p H P 1 


( 16 ) 


(17) 


(18) 


■ exp(ap H P^) If Barus' viscosity Is used 
* exp j[ln(n 0 ) ♦ 9.67][-l ♦ (1 ♦ 5 . 1 xl 0" 9 p H P 1 ) Z ]| 


If Roeland's viscosity Is used 


„ 3 2 U 

K " 4 * 1 


(19) 


( 20 ) 

(21) 


where Hq contains the term -In (8WR 2 /*)/4 Introduced In equation (3). The 
parameter Z In Eq. (20) can be expressed In terms of a and no as In 
reference 6. Furthermore, the constant load Is also Introduced by 



The unknowns Is this problem are 


( 22 ) 


*end outlet meniscus distance defined In figure 1 
N number of nodes to use as defined In figure 1 


Hq constant to define 

p e H e value of pH where dP/dX - 0 

Pj pressure at node J (J ■ 2, N) 
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The boundary conditions are 


(23) 

(24) 


p l ■ 0 for Xi - X m1n 
P - dP/dX .0 for X - X end 

Figure 1 Illustrates the calculation of X end and N, where X N Is the 
nearest node to X en( j such that X^ < X end ; H( X) can be defined as a second- 
degree curve using H^_i , H^, and Hn+i; and x end 1s calculated such that 
H(Xend) * Pe H e- Having defined X end and the remaining N ♦ 1 unknowns 
p e H e , Hq, and ?2 to P N are calculated by using an Iterative Newton- 
Raphson scheme. If the superscripts n and o are ■ -ed to define the new 
and old values of the unknowns corresponding to two successive Interatlons, we 
have 


(? e H e ) n - (? e H e )° ♦ [»(; e H e )] n (25) 

pJ ■ P° ♦ (iPj)" (26) 

H o • H o * < AH o>" (2,) 


where [a(p e H e )] n , [a(Pj)3 n , and [a(Hg)] n are now the unknowns 
to the problem. They must all be small If the convergence Is obtained. For 
the data presented, an accuracy of 1/1000 was required on each unknown. From 
the definition of the Newton-Raphson algorithm, we have for each node 1 



where (af ^/a(p e H e ) ]°, (af^/aPj) 0 , and ( af ^ /3Hq)° are defined analytically 
In appendix C. The constant load Is taken Into account by writing 

r*e nd f X end 

( aP) n dx - j - I P° dx . (AW) n 

X m1n 

(29) 

£ c, (4P,) n . <»W)" 

3-2 J J 

where Cj are the weighting factors defined In appendix 0. Since X^ 
does not coincide with X end , minor corrections are applied on the last 
values of Cj, as described In appendix B. A linear system of N *■ 1 equa- 
tions Is therefore to be solved 


J 


min 


or 
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A Gauss elimination with a partial pivoting method Is used to solve this linear 
system. 


RESULTS AND DISCUSSION 

The results of employing these techniques are shown In figures 2 to 3. 

Most of these results have been obtained by using a two-node formulatlo? for 
dP/dX because the convergence was found to be easier. 

Figure 2 gives the pressure profiles and film shapes at Iterations 0, 1, 
and 14. The 14th Iteration Is the converged solution. Barus' pressure- 
viscosity expression was used and W - 2.0452xl0 -5 , U - I.OxlfH 1 , and 
G » 5007. Barus 1 viscosity was used as this stage to enable direct comparison 
with the results of Hamrock and Jacobson (ref. 3). Starting with a Hertzian 
pressure profile (Iteration 0), the present approach converged quickly. At 
the first Iteration a pressure spike was formed that was close to the final 
converged spike. The Initial film level (Iteration 0) was obtained from 
Hamrock and Jacobson (ref. 3). In general, It took about 15 Iterations to 
obtain a converged solution, or about 2 min of CPU time on the IBM 370 com- 
puter with a mesh of 181 nodes. Contrast this with the 100 min normally taken 
by the Hamrock and Jacobson approach. At high loads Iterations have to be 
performed by starting from a previous pressure profile obtained at a lower 
load. 

Figure 3 shows the effects of a number of uniform nodes on H ra ^ n , H e , 
and the pressure spike P s . The parameters W, U, and G were exactly the 
same as those described for figure 2, and Barus 1 viscosity was also used. The 
actual pressure profiles and film shapes for four nodes are shown at the top 
of figure 3. The value of H„^ n did not change much as the number of nodes 
N max Increased: H m1n for N^x - 51 was about the same as H m ^ n for 
N max * 321 . With Increasing number of nodes H e approached an asymptotic 
value. For this case little change was observed when N^,, became larger 
than 221. When the three-node formulation of dP/dX was used, the convergence 
was more difficult to obtain, but some successful runs clearly Indicated that 
considerably fewer nodes were required for a given accuracy of H e . 
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The value of the pressure spike fluctuated, especially at low values of 
N max* because of the position of the nodes relative to the position of the 
spike. For example, at N^x ■ 51 the spike did not develop because of the 
large spacing of the nodes (sketch at top of fig. 3). Another example Is given 
by the asterisk In figure 3, for N^x ■ 241. This asterisk corresponds to the 
value of P s when a translation of half an Interval has been Imposed on the 
mesh and gives a range of fluctuation of P s . Note that an Increase of 20 was 
chosen for N^x And that more oscillations In P s might have been found If 
N max would have been varied more slowly. As the number of nodes was Increased, 
the spike became better defined. At N^x - 261 the value of P s started to 
become stable, although It was still slightly Increasing as N MX Increased. 

It can therefore be concluded that the amplitude of the spike Is limited once 
a sufficient number of nodes have been used. As the number of nodes Increased, 
there could also be large changes In the pressure spike but little correspond- 
ing change In the film thickness. The results presented In figure 3 are for a 
given uniform mesh. When a nonuniform mesh was used with many nodes near the 
spike, the pressure spike became better defined. 

Hamrock and Jacobson's (ref. 3) and the present film shapes and pressure 
profiles are compared In figure 4. Parameters were Identical to those 
described for figures 2 and 3, and N^x was fixed at 321. The value of H e 
as calculated In reference 3 Is about 2 percent larger than the value cal- 
culated In the present approach; the value of H m ^ n Is about 9 percent 
smaller. The pressure spikes also differ In amplitude and position. Although 
the differences are not yet explained, they could be attributed to the fact 
that Hamrock and Jacobson were solving the second-order differential Reynolds 
equation, whereas the present approach solved the first-order differential 
Reynolds equation. Furthermore, as noted previously, elasticity calculations 
also differ. The present results are believed to be more accurate than those 
described In reference 3, where In front of tne pressure spike region the sign 
of dP/dX was not equal to the sign of (pH - p e H e ). Furthermore, the mesh 
flow was not constant In reference 3, whereas the first-order differential 
Reynolds equation solved In the present approach Is directly representative of 
a constant flow. 

Figure 5 shows the pressure profile and the film thickness ratio H/H e 
for five dimensionless loads varying over two orders of magnitude (from 2xl0~ 5 
to 3xl0 -3 ). This dimensionless load range corresponds to a maximum Hertzian 
pressure of 0.4 to 4.8 GPa. A nonuniform mesh was always used In the spike and 
Inlet regions for the high-load cases. For those results Roelands' pressure- 
viscosity formula was used, and the dlmenslonles! speed and materials param- 
eters were fixed at U > lxlO" 11 and G « 5007. /.s the loa' 1 increased 
(fig. 5(a)), the spike became smaller and moved toward the outlet. Further- 
more, as the load Increased, the Inlet meniscus moved toward the abscissa 
X ■ -1 . The nip film thickness width and length (fig. 5(b)) both decreased as 
the load Increased. This effect could be used to predict, from experimental 
observations of the nip, the existence and amplitude of the pressure spike. 

Figure 6 shows the effect of the dimensionless load W on H e , H n « n , and 
P s . Conditions were similar to the ones defined In figure 5. At low loads 
H ffl i n was smaller than H e , but at high loads the reverse was true. The 
slopes of the curves plotted on a log-log scale are -1.16 for H e and -1.11 
for H m1n< The slope of H m ^ n is similar to the one found by Hamrock and 
Jacobson In reference 3 where a slope of -0.11 was found when plotting h m ^ n /R 
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versus W. As the load Increased, the height of the spike decreased until 
W - - 3xl0~*, where the spike no longer appeared for these given speed and 
material conditions. With a much finer nodal structure the spike might be 
observed even at higher loads, but Its height would Indeed be quite small. 

Figure 7 shows the effect on the pressure and film thickness of using 
Barus' or Roelands' pressure-viscosity expression. The dimensionless load, 
speed, and materials parameters were held fixed at W > 2.0452xl0“ 5 , 

U » lxkH 1 , and G ■ 5007. A uniform mesh of 321 nodes extending from -4 
to 1.3 was used. The large number of nodes was used to reduce the possibility 
of spike height Instabilities. Barus' formula produced a higher spike and a 
larger film thickness than Roelands 1 formula. The Increase In H e was about 
3 percent. 

Figure C shows the effect of lubricant compressibility on film thickness. 

The parameters U and G were held fixed while W was equal to 3x10“ 3 . In 
the compressible case the film thickness H e was reduced by the average density 
ratio, and H m ^ n was larger than H e . This result Is due to the effects 
of compressibility since In the Incompressible casr H m ^ n was smaller than 
H e , as wouk be expected In order to keep the mass flow constant at the nip, 
where dP/dX was large and negative. Note also that the Incompressible film 
shape was very close to the Hertzian shape (except near the nip) because of 
the high loads. 


CONCLUSIONS 

The film thicknesses and pressures In elastohydrodyna .ilcally lubricated 
contacts have been calculated for a line contact by using an Improved version 
of Okamura's approach. The new approach allows fc lubricant compressibility, 
the use of Roelands* viscosity, a very general mesh (nonconstant step), and 
very accurate calculations of the elastic deformations. The results Indicate 
that the new method Is fast and accurate. The accuracy of the film thickness, 
pressure, and pressure spike Is fully controlled by the type of mesh and the 
number of nodes used. The new approach Is not load limited as t<ere previous 
approaches. Successful runs have been made at a maximum pressure of 4.8 GPa 
with low CPU times. This approach can be a useful tool In studying rolling 
bearings and gears. Calculations of the rolling traction force and pressure 
spike at high loads are now possible and should be of first Importance for 
predicting rolling-bearing friction and fatigue life. 
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APPENDIX A 


CALCULA1 ION OF ELASTIC DEFORMATIONS 

The details of evaluating equation (6) are presented In this appendix. 

The Interval [X m ^ n , X en( j]_can be divided Into small Intervals [Xj_ lf Xj f -|] 
so that the deformation 4^4 at node 1 Is the sum of all of the small elemen- 
tary deformations dA^j calculated at node 1 and due to the pressure defined 
In the Interval [Xj_-|, Xj+-|]. 


N-l 

J. ■ 5^ dJL* (30) 

1 J-2,4,6.... 

In these small Intervals, dP/dX' Is assumed to vary linearly with X', and X* 
varies between Xji - Xj and Xj + i - Xj as Indicated on figure 9. When 
these small Intervals are used, tne distance X^ - Xj, rather than simply X^, 
has to be Introduced. 

rVr x j 

d? 1J = ‘ h J dX' (X 1 ■ x j " x,) [ ln(x 1 - X J - x ') 2 - 2 ] dX ‘ ♦ constant 

X J-1' X J 

(31) 

The linear expression for dP/dX ' reads 
dP 

^7 * (a-|X‘ ♦ a 2 ) Pj .j ♦ (a 3 X‘ ♦ a 4 ) Pj ♦ (a 5 X' ♦ a 6 ) P J+1 (32) 


where 
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a i ’ 5^ 


A 


, a . 




n 


d i ' (X j-i - V< x j-i - Vi> 


*3 ' T 2 


_ - ( Vl * Vl 1 V 

a 4 d, a 3 j f 


d 2 3 < X j - X J-1>< X J - X M> 


a 5 " < 3 


-(X. ♦ X. .) 
a, = ^ — + a r X 


(33) 


5 J 


3 = (X J*1 “ A J-l nA ;)*l 


X 4 i)(X 4 


x jV 


From the expression for dP/dX', It can be seen that can be expressed as 


d *1j * d0 1,J-l P J-1 * dD 1,J P J * dD 1.J*l P J*l + constant 
where dD^j are the elementary Influence coefficients calculated as 
.X, ,-X, 


(34) 


dD 


p»l 1 

i.j-1 * - hj <V * v« x , - x j - x '>[ 1n < x i - x 3 - x ') 2 - 2 J 


dX' 


X M' X J 


(35) 


A similar relation that uses the corresponding coefficients 83 , a+ and 
35 , 35 Is used to define dD^ i and dD^ hi . By adopting the change of 
variables 


Z 3 X 1 - X J - X *' Z m1n 3 X 1 - Vl ; Z max 3 X 1 ' X j*l 

b 2 = a 1 (X 1 - X^) ♦ a 2 ; dZ = - dX' 


(36) 


dDi j_i Is calculated as 
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dD i.M ■- b ,*♦■*,) z [in <z 2 ) - 2 ] « 
Z m1n 

/ Z max /* Z max 5 

Z In Z dZ ♦ a 1 / Z* In 

*m1 n Z m1 n 


Z‘ dZ ♦ 2b, 



z 


Z 1 

z 2 

2 

max 

- 2a l 

z 3 

3 

max 


Z m1n 


Z m1n _ 


* • U 


Z^ 2 7 ^ / 3 \ Z Z Z^ 

-2b 2 J- (In z' - 1) ♦ 2a^ ^^ln |Z| J - l) ♦ 2b ? f- - 2a, f- 


z max 

Z m1n 

(37) 


Using the variables M, N, and K, we obtain after rearrangement 

dD, 


'1 

dD 


.J-l ’ ‘ 2* ( 
1.3 * ‘ ii (‘ 


a-jK ♦ a 2 


**- ' a 3 K + a 4 2 


(38) 


dD 


1J*1 * ' 2* ( a 5 K f a 6 2 ) 


where 


" ■ l L ( ,n Z L - 3 ) - 4, ( ,n 4« - 3 ) 


N - Z 


max 


( ln 'W 3 - 4 ) - 4n( ,n 


'W - 4 


) 


(39) 


There Is no problem concerning the singularity occurring for ■ Xj. 

When there Is a singularity In the Interval X^-|] ( we can Integrate 

along the two half Intervals [X^_-|, X^] and [X^, X^-j ] and use the relation 

11m (Z 2 ln Z 2 ) » 11m (Z 3 ln |Z| 3 ) - 0 (40) 

Z-»0 Z-»0 
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which shows the relations [W) to be also valid for ■ Xj. Finally, the 
deformation 4^ Is obtained by summing all of the elementary small deforma- 
tions di^j 


N-l 

£ 

J-2,4,6, 


di 


1.J 


l 


1 




Y/j 


- 4 1n 



(41) 

(42) 


where D^j - dO^ < If J Is even. If J Is odd (J ■ 3 In the following 
example), 'we calculate dD ^3 coming from the Interval [X], X 3 ] and add to It 
the value dD^j coming from the Interval [X 3 , X 5 ] In order to obtain the final 
value of 0 ^ 3 . For the first and last values of D^j, we simply have 


° 1,1 * d ° 1 .1 


(43) 


°1,N a d °1,N 

Note now that D. *Pi Independent of the load. At high loads the 
j ' »J J 

dimensionless film thickness H e becomes very small with respect to 4^. 
However, by us^lng an appropriate change of variable (X 3 x/bc), the maximum 
deformation 4 ma - can be kept equal to H e . Using the last change of vari- 
able, we have 


*, ■ c £ Y p j - « ,n ( r2 t* c? ) < 44 > 

where are the new Influence coefficients obtained with the new value 

of X. From the definition of c we have 

<«> 

The maximum deformation 4 max Is close to the maximum Hertzian deformation 4 h 
( defined In eq. ( 10 )), leading to 

(r^c 2 ) (46) 

The value of c can therefore be defined as 

c = ^ exp (2H e - 0.5) (47) 

This numerical "trick* was definitively helpful at high loads but was not used 
In the results presented herein. 
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APPENDIX B 


CORRECTIONS TO BE APPLIED TO WEIGHTING FACTORS DUE TO X 

For X - X en( j , P - dP/dX • 0. It Is assumed that between X^ and X en( j, 
p( X) Is described by a second-degree polynomial. To respect the previously 
mentioned boundary condition, we have 


< x - w 

" ( X N - x end) 


2 N 


d£ 

dX 


For X » X|| we have therefore 


2(X - X_J 


end 

( X N * X end) 


- P 
2 r N 


(♦ 8 ) 


(49) 


X N - X end N 


(50) 


which can also be expressed as 

Vn-l ■ 0: ‘ 


n,n 


x - X 
*N end 


a . ■ 0 
n,n+l 


The Integral of the pressure between X^ and X en( j leads to 

•X 


/ 


end 


p dx * t l*end - V P N 


aC N P N 


(51) 


(52) 


If N Is odd, we substract from the value of Cf| the value dCfj+i,N coming 
from the Interval [X|j,Xn+ 2] and add aC N to set the final value of* Cr. If 
N Is even, we modify C^.i and C N . From C N _i , we substract dC|| v n_i 
from the Interval (Xjj_-|, Xn+i ] and add the weighting factor ACr_i coming from 
the Integration of P between Xr_i and X^. Using the "trapeze" rule, 


AC 


X - X 

a n *N-1 


N-l 


The value of C|j Is finally defined as 

C N " &C N-1 


♦ AC. 


(53) 


(54) 
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APPENDIX C 


CALCULATION OF JACOBIAN FACTORS 
The factors to be defined are 

afi 


1 


af l af l 

»(?, H e ) : ’Y W 0 


with 


f h 3 t&\ k; I H 

f 1 " H 1 Idxj . ‘ Kn 1 l H 1 " - 

\ / 1 \ p. 


(16) 


H 


1 " H 0 * 2 1 * * ^ °1J P J 


(17) 


dP/dX, and p^ are defined by equations (11), (19) or (20), and (18) 
Before the final Jacobian factors are defined, we can define 


3H. 3H 

8P^ ' °1J 


A 


3n 


1 


apj " ° P H n 1 k 1J 


3n 


If Barus' viscosity Is used 

Z-l 


M <J _Q - 

= 5.1x10 * P H (In n 0 ♦ 9.67)(1 ♦ 5.1x10 ’p^) 

J 

If Roelands' viscosity Is used 

where k^ Is the Kronecker symbol (k^j - 1 If 1 ■ J; « 0 If 1 * J) 


3(1/^) 


0.6xl0“ 9 p 


H 


3P 


J 


1 ♦ 2.3xl0" 9 p H P 1 


3Pi 




“ a 1,1-l k 1-l,J * a 11 k 1J * a 1,Hl k Hl,J 


(55) 


J 


It Is now easy to define the Jacobian factors 
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af 1 Kn, 

a(p e H e ) Py 


af 

ap 


1 

J 



af 

3H 


1 

o 





(56) 
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APPENDIX D 


DEFINITION OF WEIGHTING FACTORS Cj 

The weighting factors Cj are defined by using the three-point Lagrange 
polynomial with a general mesn (nonconstant step). In the Interval [Xj_-| , 

Xj f i } the pressure Is described by a second-degree polynomial In X 1 

. x,(x ' p . ii : x i - x i-i> <«■ * s - x u!> p 

<x« - X. ,)(X. - X.,.) r i 


(*M - *,)(*!.! - x J*i> J -’ 


<*' * x i - x i-i> x ' 


<x 3»i * x j-l )(x j*l ' V 

We can now define the coefficients dc j,j-i» dCjj, and dCjj f -| such * hat 


(57) 


/ 


X' 

max 


P dX ‘ * dC J,J-l P J-l + dC JJ P J * dC J,J*lVl 


min 


dC 


j.j-i - x j )( Vi - Vl> 


,3 


dC 


J,J (Xj - Xj_ 1 )(Xj - x Jfl ) 


,3 


;• 

max 

< X M 

,xi 2 
' 2 

X ™x 

X' 

min 


-.2 

x j.l> 

X' 

2 


\ 


(59) 


* (X J ‘ x j_ 1 )(Xj - x J+1 ) X' 


X' 

max 


X' 

A m1n 


dC 


J.J* 1 (Xj +1 - x j_i ) ( " x j> 


w 1 3 x ,2 

r ♦ (X j - x j-i> * 


X' 

max 


X' 

min 


) 


These coefficients are calculated for J » 2,4,..., N - 1. The coefficients 
Cj are finally defined as 


Cj-dC JtJ If j is 


even 
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When J Is odd (J - 3 In the following example), we calculate dC? 3 corre- 
sponding to the Interval [X3 v X3 ] and add to It dC^ coming from tfie Interval 
[X3,Xc] In order to define Cj or C3 In this example. Minor corrections 
are also applied on the last values of Cj In order to respect the boundary 
condition for X ■ X en( j, as shown In appendix B. 
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Fig. 1 Sketch to illustrate calculations of X end and N. 


PRESSURE. P; AND FILM SHAPE. H/2H, 



X-COORDINATE, x/b 


Fig. 2 Pressure profiles and film shapes at iterations 0. 
1, and 14. Barus 1 pressure-viscosity formula; W ■ 
2.0452x10-5; U -l.OxlO-H; and G - 5007. 
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Fig. 3 Effect of number of uniform nodes on H n s n , H e , and 
P. for a fixed W ■ 0. 20452xl0" 4 , U - lxlO -11 ,’ dnd G ■ 5007. 


Bams 1 pressure-viscosity formula assumed. 
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Fig. 4 Comparison between Hamrock and Jacobson's pres-* 
sure profile and film shape (3) and the authors'. Barus' 
pressure-viscosity formula; W • 2. 0452xl(T* ; U - 
l.OxlO’ 11 ; G -5007; and N mav - 321. 



FILM SHAPE. H/H e 



H) Pressure profiles. 



(b) Film shapes. 

Fig. 5 Pressure profiles and film shapes for varying 
dimensionless loads. Roelands 1 pressure-viscosity 
formula; U-1.0xl0" 11 ; and G- 5007. 



DIMENSIONLESS FILM THICKNESS, 



DIMENSIONLESS LOAD, W 

Fig. 6 Effect of dimensionless load W on H^, H m j n , and 
° s . Roelands' formula for pressure-viscosity assumed; 
U -l.OxlO" 11 ; and G -5007. 



PRESSURE. P; OR FILM SHAPE. H/2(He) R 



Fig. 7 Effect of using Barus' or Roelands' pressure- 
viscosity formula on pressure and film shape. Dimen- 
sionless load, speed, and materials parameters held 
fixed at W • 0. 2045xl0" 4 . U - lxlO' 11 , and G - 5007. 
Uniform mesh of 321 modes extending from X equal 
-4 to 1.3 used. 



incompressible 


1.4 


incompressible 

compressible 



Fig. 8 Effect of lubricant compressibility on film thickness. 
U » 1.0xl0 _ ll; W - 3x10-3; and G ■ 5007. 



i.j+1 



Fig. 9 Calculation of deformation d6j; at node i due to pressure 
acting in interval IXj.j, Xj + 1 l. J 
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